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ABSTRACT 

A  new  approach  is  developed  to  reconstruct  a  three-dimensional  incompressible  flow  from  noisy  data  in  an 
open  domain  using  a  two-scalar  (toroidal  and  poloidal)  spectral  representation.  The  results  are  presented  in  two 
parts:  theory  (first  part)  and  application  (second  part).  In  Part  I,  this  approach  includes  (a)  a  boundary  extension 
method  to  determine  normal  and  tangential  velocities  at  an  open  boundary,  (b)  establishment  of  homogeneous 
open  boundary  conditions  for  the  two  potentials  with  a  spatially  varying  coefficient  k,  (c)  spectral  expansion 
of  k,  (d)  calculation  of  basis  functions  for  each  of  the  scalar  potentials,  and  (e)  determination  of  coefficients  in 
the  spectral  decomposition  of  both  velocity  and  k  using  linear  or  nonlinear  regressions.  The  basis  functions  are 
the  eigenfunctions  of  the  Laplacian  operator  with  homogeneous  mixed  boundary  conditions  and  depend  upon 
the  spatially  varying  parameter  k  at  the  open  boundary.  A  cost  function  used  for  poor  data  statistics  is  introduced 
to  determine  the  optimal  number  of  basis  functions.  An  optimization  scheme  with  iteration  and  regularization 
is  proposed  to  obtain  unique  and  stable  solutions.  In  Part  II,  the  capability  of  the  method  is  demonstrated  through 
reconstructing  a  2D  wind-driven  circulation  in  a  rotating  channel,  a  baroclinic  circulation  in  the  eastern  Black 
Sea,  and  a  large-scale  surface  circulation  in  the  Southern  Ocean. 


1.  Introduction 

Ocean  observational  current  data  are  usually  acquired 
from  limited  number  of  stations  in  domains  with  open 
boundaries  and  contain  various  errors  or  noises.  It  is  an 
important  task  for  physical  oceanographers  to  establish 
(or  to  reconstruct)  a  realistic  and  complete  velocity  field 
from  sparse  and  noisy  data. 

From  a  mathematical  point  of  view,  the  reconstruction 
requires  solving  a  least  square  problem  without  or  with 
a  priori  information  (limit)  on  the  circulation  charac¬ 
teristics.  An  a  priori  limit  can  be  formulated  as  a  set  of 
inequalities  that  the  solutions  should  satisfy,  as  a  dy¬ 
namical  model  applied  to  the  description  of  circulation 
dynamics  or  hypotheses  on  statistical  properties  of  re¬ 
constructed  field. 

Several  techniques  are  available  for  fulfilling  such  a 
task:  various  kinds  of  spline  interpolation  (e.g.,  Washba 
and  Wendelberger  1980;  Smith  and  Wessel  1990;  Bran- 
kart  and  Brasseur  1996),  optimal  interpolation  (OI;  e.g.. 
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Gandin  1965),  fitting  models  (Cho  et  al.  1998;  Lipphardt 
et  al.  1977,  2000),  objective  mapping  combined  with  a 
fitting  (e.g.,  Davis  1985),  and  numerous  approaches  us¬ 
ing  ocean  numerical  models  such  as  the  adjoint  method 
and  Kalman  filter  (e.g.,  Malonette-Rizzoli  and  Tziper- 
man  1996). 

Several  error  sources  deteriorate  the  reconstruction 
skill.  One  of  them  is  the  uncertainty  in  boundary  con¬ 
ditions,  especially  at  open  boundaries.  Therefore,  how 
to  determine  open  boundary  conditions  becomes  a  key 
issue  in  the  reconstruction  process. 

The  classical  OI  technique  does  not  allow  accounting 
for  any  boundary  condition  as  an  a  priori  limitation. 
Davis  (1985)  used  a  combined  Ol-spectral  fitting  model 
with  a  priori  knowledge  of  the  statistical  weights  to 
overcome  this  weakness.  However,  it  remains  uncertain 
how  to  select  the  weights  for  an  open  domain  and  how 
to  determine  basis  functions  with  a  priori  nonzero  flux 
at  the  open  boundary. 

With  velocities  given  along  the  open  boundary  and 
with  an  additional  boundary  condition  such  as  the  “nat¬ 
ural”  boundary  condition  (Courant  and  Hilbert  1966), 
the  spline  functions  can  be  used  as  universal  basis  func¬ 
tions.  However,  a  detailed  analysis  (Inoue  1986)  shows 
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that  the  natural  boundary  condition  is  more  appropriate 
for  rigid  than  open  boundaries. 

An  effective  approach  to  determine  open  boundary 
condition  is  to  use  an  optimization  procedure  together 
with  an  ocean  numerical  model  such  as  the  adjoint  meth¬ 
od  (e.g.,  Seiler  1993),  the  Jacobian  matrix  method  (e.g., 
Chu  et  al.  1997),  and  the  local  variation  principle  (e.g., 
Shulman  and  Lewis  1995).  The  open  boundary  condi¬ 
tions  can  be  determined  through  minimizing  the  differ¬ 
ence  between  observations  and  model  estimations  inside 
the  integration  domain.  These  approaches  need  to  solve 
an  ill-posed  inverse  problem  using  an  iteration  proce¬ 
dure  independent  of  the  chosen  local  boundary  condi¬ 
tions.  To  ensure  their  convergence  to  physical  reality, 
we  need  a  “good”  initial  approximation  for  the  open 
boundary  conditions  even  if  the  noise-to-signal  ratio  in 
the  data  is  small  (Engl  et  al.  1996). 

Without  knowing  statistical  weights  and  without  us¬ 
ing  ocean  numerical  models,  a  kinematical  method  is 
proposed  for  reconstructing  a  velocity  field  from  noisy 
and  sparse  data.  For  a  three-dimensional  incompressible 
flow,  two  scalar  functions,  toroidal  ('VT)  and  poloidal  (d>) 
potentials,  satisfy  Poisson  equations  with  the  vertical 
vorticity  (£  =  dvldx  —  du/dy )  and  vertical  velocity  (w) 
as  the  sources  terms,  respectively  (Moffat  1978;  Ere- 
meev  et  al.  1992a, b;  Chu  1999;  Chu  et  al.  2002). 

The  two  potentials  (<E>,  'T)  are  expanded  into  a  series 
of  basis  functions,  which  are  eigenfunctions  of  the  La- 
placian  operator  for  the  given  geometry  of  the  ocean 
domain  and  homogeneous  boundary  conditions.  For  a 
closed  basin  with  no  slip  on  the  rigid  boundary,  it  is 
typical  to  use  the  Dirichlet  and  Neumann  boundary  con¬ 
ditions  for  and  $  (Eremeev  et  al.  1992a,b),  and  to 
use  the  Dirichlet  boundary  condition  for  the  transport 
streamfunction  (Rao  and  Schwab  1981).  This  procedure 
allowed  construction  of  basis  functions  and  reconstruc¬ 
tion  of  the  velocity  field  from  data  for  the  Black  Sea 
(Eremeev  et  al.  1992a, b)  and  Lake  Ontario  (Rao  and 
Schwab  1981),  both  treated  as  closed  basins. 

Lipphardt  et  al.  (2000)  proposed  an  approach  for  a 
domain  with  an  open  boundary  through  expanding  the 
poloidal  potential  into  a  combination  of  harmonic  func¬ 
tions  satisfying  boundary  conditions  a  priori  given  from 
a  large-scale  numerical  model  and  a  spectral  decom¬ 
position  with  the  basis  functions  as  the  eigenfunctions 
of  the  Laplacian  operator  with  homogeneous  Neumann 
conditions.  The  reconstruction  skill  of  their  method 
strongly  depends  on  the  quality  of  the  numerical  model. 

In  this  study,  a  new  set  of  basis  functions  is  introduced 
for  reconstructing  the  ocean  circulation  in  a  domain  with 
open  boundaries.  These  functions  are  the  eigenfunctions 
of  the  Laplacian  operator  with  homogeneous  mixed  con¬ 
ditions.  With  known  velocities  along  the  open  boundary, 
the  mixed  boundary  conditions  are  accurate.  With  un¬ 
known  velocities  along  the  open  boundary,  a  parame¬ 
terization  scheme  is  proposed  for  obtaining  approximate 
open  boundary  conditions  from  data.  In  general,  the 
reconstruction  is  reduced  to  linear  and  nonlinear  re¬ 


gression  models  for  known  and  unknown  velocities 
along  the  open  boundary,  respectively.  For  the  latter 
(without  data  on  the  open  boundary),  the  velocity  inside 
the  domain  and  along  the  boundaries  are  simultaneously 
determined. 

The  reconstruction  skill  depends  on  various  factors 
such  as  the  noise  level  in  the  data,  the  sampling  strategy, 
and  the  quality  of  numerical  algorithms  used  in  the  re¬ 
construction.  This  study  shows  that  even  if  the  com¬ 
bination  of  all  the  factors  is  not  favorable,  the  recon¬ 
structed  circulation  pattern  and  the  velocity  along  the 
open  boundary  are  acceptable  according  the  criterion 
mentioned  in  the  text.  Thus,  the  proposed  reconstruction 
scheme  is  a  useful  tool  in  analyzing  kinematic  and  en¬ 
ergetic  characteristics  of  the  circulation,  and  provides 
zero-order  initial  conditions  for  numerical  ocean  mod¬ 
els. 

The  outline  of  this  paper  is  as  follows.  In  section  2, 
we  describe  the  decomposition  of  flow  into  poloidal  and 
toroidal  components,  specify  the  boundary  conditions 
for  both  potentials,  and  introduce  two  sets  of  complete 
basis  functions.  In  section  3,  we  depict  the  principal 
relationships  used  in  the  reconstruction  process.  In  sec¬ 
tion  4,  we  discuss  a  regularization  technique.  In  section 
5,  we  provide  an  approach  to  determine  the  optimal 
number  of  basis  functions  used  in  linear  or  nonlinear 
regression  models,  and  two  criteria  such  as  Vapnik- 
Chervonkis  cost  function,  and  root-mean-square  errors 
to  determine  the  reconstruction  skill.  In  section  6,  we 
give  conclusions. 

2.  Two  scalar  potentials 

a.  Toroidal  and  poloidal  components 

In  magnetohydrodynamics  and  astrophysics,  it  is 
common  to  decompose  any  vector  Q  in  an  arbitrary 
coordinate  system  into  three  parts  (Dubrovin  et  al. 
1992).  For  example,  it  is  in  spherical  coordinate  system, 
written  as 

Q  =  r  X  VA1  +  rA,  +  VA3,  (2.1) 

where  A,,  A2,  and  A3  are  scalar  functions  (see  appendix 
A),  V  is  the  three-dimensional  gradient  operator,  and  r 
is  the  radius  vector  from  the  origin.  Borrowing  this  idea 
for  ocean  currents  (Q  =  u)  satisfying  the  incompressible 
property 

V  u  =  0,  (2.2) 

the  three-dimensional  velocity  field  at  large-,  meso-,  and 
submesoscales  is  represented  by  [(A. 5),  see  appendix 
A]: 

u  =  V  X  (rTO  +  V  X  V  X  (rd>),  (2.3) 

where  the  two  terms  in  the  right-hand  side  of  (2.3)  are 
called  toroidal  and  poloidal  velocities. 

If  the  velocity  is  reconstructed  on  horizontal  planes, 
the  radius  vector  r  can  be  replaced  by  the  unit  vector 
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Fig.  1.  Boundary  conditions  and  open  boundary  extension  method. 


where  (n,  r)  represent  normal  and  tangential  unit  vectors 

along  r  u  r;. 


b.  Rigid  boundary  segment  ( T) 

Along  the  rigid  boundary  ( T )  segment,  the  kinematic 
boundary  conditions  are  given  by 


V„  =  0,  PT  ^  0. 


The  two  scalar  potentials  satisfy  (Eremeev  et  al. 
1992a,b) 


T'lr  =  C, 


d2<t> 

dndz 


=  0, 


(2.7) 


in  the  vertical  direction  k  (Moffat  1978).  Thus,  the  ve¬ 
locity  u  ( u ,  v,  vv),  determined  on  any  horizontal  plane, 
is  represented  by  (Eremeev  et  al.  1992a,b) 

ibV  d'V  d23> 

u  = - 1 - ,  v  = - 1 - ,  and 

By  dxdz  dx  dydz 

w  =  -  Ad>,  (2.4) 

where  the  Cartesian  coordinate  system  is  used  with  (x, 
y)  and  z  as  the  horizontal  and  vertical  coordinates,  re¬ 
spectively. 

Obviously,  both  toroidal  and  poloidal  potentials  sat¬ 
isfy  the  Poisson  equations 


where  C  is  a  constant  to  be  determined.  Since  ra  and 
rb  are  also  the  two  end  points  of  the  rigid  segment,  we 
have 

'i'  I  Ta  =  ^  I  Tt  =  C. 

Note  that  the  scalar  potentials  "T  and  <b  cannot  be  de¬ 
termined  if  the  velocity  at  the  rigid  segment  vanishes 
(Ladyzhenskaya  1969): 

K  =  vT  =  o. 


c.  Open  boundary  segment  ( f]) 


A^f  =  -£  and  (2.5a) 

=  -w.  (2.5b) 

Here,  A  =  d2/dx2  +  d2/dy 2  is  the  two-dimensional  La- 
placian  operator,  and  £  =  dv/dx  —  du/dy  is  the  vertical 
component  of  vorticity. 

In  general,  the  toroidal  and  poloidal  potentials  Jp  and 
d>  are  not  the  same  as  the  geostrophic  stream  function 
and  velocity  potential  commonly  used  in  meteorology 
and  oceanography  (e.g.,  Lynch  1988).  If  the  Coriolis 
parameter  varies  considerably  within  the  domain,  the 
poloidal  potential  satisfies  the  Poisson  equation  with  a 
source  term  determined  by  the  horizontal  velocity  and 
the  gradient  of  the  Coriolis  parameter  even  in  the  pure 
geostrophic  flow.  That  can  be  checked  out  through  the 
direct  substitution  of  (2.4)  into  the  geostrophic  equa¬ 
tions. 

Before  introducing  a  new  set  of  basis  functions  for 
an  open  domain,  one  should  answer  the  question,  What 
boundary  conditions  will  be  used  for  toroidal  and  po¬ 
loidal  potentials?  Let  fluid  be  in  a  simply  connected 
domain  with  rigid  (T)  and  open  (T[)  boundaries.  The 
open  boundary  has  two  end  points  r„  and  rh  (Lig.  1). 
Both  toroidal  and  poloidal  potentials  satisfy  Poisson 
equations  with  coupled  boundary  conditions.  Along  the 
boundary,  the  velocity  is  decomposed  into  normal  (V„) 
and  tangential  ( VT )  components 

d2d>  rTP  d2<J) 

V„  =  —  +  - ,  Vr= - +  - ,  (2.6) 

dr  dndz  dn  drdz 


When  the  normal  ( V„ )  and  tangential  (VC)  components 
are  given,  we  integrate  the  first  equation  in  (2.6)  with 
respect  to  r  along  the  open  boundary  segment  (T[)  from 
the  end  point 


'P(t)  =  )p|Ti 


d2d>\  , 

- I  dr  , 

dndz) 


(2.8) 


and  define  a  coefficient  k(t)  varying  along  TJ  by 


k(t) 


dr’  +  C 


(2.9) 


The  boundary  condition  for  ''P  is  obtained  from  (2.6): 


d'if 

—  +  /c(r))p 
dn 


=  0, 


r; 


(2.10) 


where  r  s  [ra,  t6],  is  a  moving  point  along  T[. 

When  the  constant  C  in  (2.9)  vanishes,  the  coefficient 
k(t)  tends  to  infinity  as  r  tends  to  ra  or  rb.  Prom  the 
theoretical  point  of  view,  such  a  behavior  of  k(t)  does 
not  add  any  complexity  in  applying  (2.10).  Pirst,  the 
singularity  of  k(t)  can  be  effectively  eliminated  by  a 
perturbation  technique  (Morse  and  Preshbach  1953)  due 
to  the  toroidal  potential 'T  being  a  smooth  function  along 
the  boundary  TJ.  Second,  we  may  transform  (2.10)  into 
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dV 

Ki(T)—  + 

on 


=  0, 


where  no  singularity  occurs.  Here 


KlO) 


dndzj 


+  C, 


k2(t) 


In  this  paper,  these  conditions  are  not  used  because 
the  unknown  open  boundary  condition  is  assumed  a 
priori.  A  smooth  function  is  used  to  parameterize  k( t) 
and  in  turn  to  provide  an  approximate  open  boundary 
condition. 


Thus,  the  boundary  conditions  for  the  rigid  and  open 
segments  are  represented  by  (2.7)  and  (2.10),  (2.13), 
respectively.  The  corresponding  basis  functions  are  the 
eigenfunctions  of  the  following  spectral  problems: 

A%  =  —A  A,  (2.14) 

A<D,„  =  (2.15) 

—  =0,  and  (2.16) 

dn  r 

<DJri  =  0,  (2.17) 


**lr  =  0, 


iFVk 

—  +  <r)9k 

on 


=  0, 


which  are  the  mixed  boundary  conditions  formulated 
for  a  simply  connected  domains. 


d.  Basis  functions 


(ii)  Features  of  the  basis  functions 


1)  Simply  connected  enclosed  domain 

For  a  domain  with  an  enclosed  rigid  boundary  (T), 
the  appropriate  basis  functions  { 'vtrit }  and  { }  are  given 
by 

k  =  1,  .  .  .  ,  oo, 

(2.11) 

m  =  1,  .  .  . ,  oo. 

(2.12) 

Here,  {At}  and  {/r,„}  are  the  corresponding  eigenvalues. 
Similar  basis  functions  were  obtained  by  Rao  and 
Schwab  (1981)  for  Lake  Ontario,  and  Eremeev  et  al. 
(1992a, b)  for  the  Black  Sea,  respectively. 

2)  Simply  connected  open  domain 

Both  velocity  components  (V„,  VT)  are  usually  non¬ 
zero  at  the  open  boundary  TJ.  The  toroidal  potential 
satisfies  the  boundary  condition  (2.10)  with  the  param¬ 
eter  k(t)  depending  on  <t>.  Thus,  a  possible  approach  is 
to  specify  at  the  open  boundary  to  obtain  a  condition 
for  T'. 

(i)  Minimum  poloidal  kinetic  energy  assumption 

Decomposition  of  the  velocity  field  into  toroidal  and 
poloidal  parts  can  take  various  forms.  The  potentials 
have  no  physical  significance  themselves  (Ladyzhen¬ 
skaya  1969).  They  are  meaningful  only  in  representing 
the  circulation.  To  reduce  the  degree  of  freedom  without 
loss  of  any  generality,  the  poloidal  kinetic  energy  is 
assumed  averaged  over  the  domain  including  the  open 
boundary  to  be  minimal  and  obtain,  according  to  Ped¬ 
ersen  (1971), 

d$ 

—  =  0  at  r;.  (2.13) 

dz 


=  -AA,  ^*lr  =  0, 


d<I> 

=  -UL  <b  — - 

m  r^m  nn  ~ 

rill 


=  o. 


The  basis  functions  defined  by  (2.14)-(2.17)  have  the 

following  features. 

1)  Each  of  the  two  sets  of  basis  functions  {''P{ }  and 
{<&„,}  is  orthonormal  and  complete  (Vladimirov 
1971).  To  calculate  directly  these  basis  functions,  it 
requires  a  priori  knowledge  of  geometry  and  velocity 
components  at  the  boundary  (i.e.,  a  known  boundary 
condition).  For  unknown  boundary  conditions,  a 
nonlinear  regression  scheme  should  be  developed. 

2)  Three  reasons  make  the  basis  functions  defined  here 
more  appropriate  than  trigonometric  polynomials 
(plane  geometry)  and  spherical  harmonics  (spherical 
geometry)  in  flow  reconstruction  from  noisy  and 
sparse  data.  First,  the  trigonometric  polynomials  and 
spherical  harmonics  are  not  the  solutions  of  (2.14)- 
(2.17)  for  a  domain  with  complex  boundaries  and / 
or  with  k  varying  along  the  open  boundary.  That  is 
to  say  that  the  trigonometric  polynomials  and  spher¬ 
ical  harmonics  cannot  formulate  a  complete  set  of 
basis  functions  in  this  case.  Second,  the  spectral  se¬ 
ries  usually  converges  quicker  using  the  basis  func¬ 
tions  determined  by  (2.14)-(2.17)  than  using  trigo¬ 
nometric  polynomials  and  spherical  harmonics  since 
the  physical  information  at  the  boundary  is  suffi¬ 
ciently  used.  This  leads  to  fewer  modes  needed  using 
{'vLlt}  and  {<&,„}  as  the  basis  functions  than  using  the 
trigonometric  polynomials  and  spherical  harmonics. 

3)  If  normal  and  tangential  velocities  along  the  open 
boundary  change  with  time,  the  coefficient  k  also 
depends  on  time.  The  velocity  field  should  be  re¬ 
constructed  at  a  particular  time.  This  usually  does 
not  add  any  complexity  to  the  reconstruction.  How¬ 
ever,  it  may  consume  considerable  computer  re¬ 
sources  if  (2.9)  is  applied  for  the  analysis  of  long¬ 
term  observation  series.  This  topic  will  not  be  dis¬ 
cussed  further  in  this  paper. 

4)  The  boundary  condition  (2.10)  may  be  simplified, 
for  example,  assuming  zero  tangential  velocity  at  an 
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extended  open  boundary  T,.  The  boundary  condi¬ 
tions  become  (Fig.  1) 


<TT 

dn 


=  0, 


5$ 

dz 


=  0. 


(2.18) 


The  solutions  of  (2.14)  and  (2.15)  are  relatively  easy 
to  obtain  since  {T^}  and  {<Fm}  do  not  depend  on 
each  other.  Notice  that  the  condition  (2.18)  allows 
for  nonzero  normal  velocity  along  the  open  bound¬ 
ary. 

5)  The  approach  can  be  extended  to  a  multiply  con¬ 
nected  domain  through  the  methodology  originally 
described  by  Kamenkovich  (1961).  This  is  demon¬ 
strated  through  circulation  reconstruction  in  the 
Southern  Ocean  from  the  drifter  data  [see  Part  II, 
(Chu  et  al.  2003)]. 


3)  Multiply  connected  open  domain 

Consider  a  multiply  connected  domain  with  a  rigid 
boundary  T  and  an  open  boundary  P.  Inside  the  do¬ 
main,  there  are  L  islands  with  rigid  boundaries  ri5  T, 
.  .  . ,  Tl  (Fig.  2).  The  solution  of  the  Poisson  equation 
(2.5a)  can  be  decomposed  into  (Vladimirov  1971): 


constants  along  the  boundaries  ri5  T,,  .  .  .  ,  F,,  and  T 

u  r\ 

The  variable  'T  satisfies 


AT'  =  (2.22) 


with  the  homogeneous  boundary  conditions, 


T^r1,r2,...,r1,r  0’ 


d'if  ~ 

-  +  /AP 

dn 


where  the  coefficient  k  is  given  by 


0.  (2.23) 


K 


d2<t> 

- C — 

drdz  dn 


(2.24) 


The  ( L  +  1)  constants,  C,(Z  =  1,  2,  .  .  .  ,  L)  and  C 
should  be  determined  by  the  data  or  by  additional  dy¬ 
namic  constraints  (Kamenkovich  1961;  McWilliams 
1977;  Flierl  1977). 

The  poloidal  potential  satisfies 


A<I>  =  —w, 


with  the  rigid  boundary  conditions 


dzdn 


=  0  on  T,,  .  .  .  ,  TL  and  T,  (2.25) 


and  with  the  open  boundary  conditions 

dT> 

—  =  0  on  T'. 

dz 


(2.26) 


The  capability  for  the  multiply  connected  domain  is 
demonstrated  through  reconstructing  the  Southern 
Ocean  circulation  from  the  First  Global  Atmospheric 
Research  Program  (GARP)  Global  Experiment  (FGGE) 
drifter  buoy  data,  which  is  illustrated  in  Part  II  (Chu  et 
al.  2003)  of  this  paper. 


3.  Reconstruction  procedure 


T'  =  2  C,/Tf(,)  +  CT'  +  T>,  (2.19) 

1=1 

such  that  the  harmonic  functions  'I'"1  and  'T  satisfy  the 
Laplacian  equation 

AT'(')  =  0,  AT'  =  0,  (2.20) 

with  the  inhomogeneous  boundary  condition 

~  f 1  on  r, 

|^0  on  r,  r\  and  Tr(l'  ^  Z), 


a.  Horizontal  velocity 

From  the  mathematical  point  of  view,  such  a  recon¬ 
struction  can  be  reduced  to  a  linear  regression  if  normal 
and  tangential  velocities  are  given  at  the  open  boundary, 
and  to  a  nonlinear  regression  otherwise.  The  regression 
coefficients  are  calculated  from  observed  data. 

The  two  scalar  potentials  are  expanded  into  following 
Fourier  series: 

T'Xx,  y,  z,  f)  =  2  ak{z,  t°)'Vk{x,  y,  z,  k°), 

*=1 


fi  on  r  u  r 

jo  on  r;  (Z  =  1,  2,  .  .  .  ,  L). 


(2.21) 


d<I>(x,  y,  z,  t°) 
dz 


2  b,„(z,  f°)T>,„(x,  y,  z), 


(3.1) 


Here,  C;(Z  =  1 ,  2,  .  .  .  ,  L)  and  C  are  the  integration  where  t°  is  the  time  of  observations. 
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Substituting  (3.1)  into  (2.4)  and  truncating  modes  at 
K  and  M  for  ^  and  <f>,  respectively,  we  obtain  an  ap¬ 
propriate  regression  model  through  minimizing  the  fol¬ 
lowing  functional 

/(«i,  .  .  . ,  aK,  bu  .  .  .  ,  bM,  k,  P) 

=  ^(ll<s  -  MjJIp  +  IKS  -  VKM\\2P)  min,  (3.2) 


where  w°bs,  u°bs  are  observations  at  location  p  (p  =  1, 
.  .  .  ,  P),  P  is  the  total  number  of  observations,  and 


UKM 


K 


2  o 

*=i 


3^t(x,  y,  z,  k°) 
dy 


M 


+  2  bmiz,  t°) 


ddbnCr,  y,  z) 
dx 


v 


KM 


K 


-2  f°) 


y,  z,  O 
dx 


M 


+  2  t° ) 


y,  z) 
dy 


(3.3) 


are  reconstructed  velocity  components.  Here,  ||.  .  ,||P  is 
the  P-dimensional  Euclidean  norm.  After  obtaining  the 
spectral  coefficients  in  (3.1)  and  parameterized  k,  the 
horizontal  velocity  can  be  calculated  for  any  point  of 
the  domain. 


b.  Vertical  velocity  and  vorticity 

For  the  reconstruction  of  vertical  velocity  at  depth 
z*,  vertical  differentiation  is  used  on  the  second  equation 
in  (2.5): 


dw  ,  dd> 
—  =  -A — . 
dz  dz 


(3.4) 


Equation  (3.4)  is  integrated  from  the  surface  to  the  depth 
z*  or  from  the  depth  z*  to  bottom  and  rewrite  the  result 
of  integration  as  the  following  integral  equation: 


G(x,  x\  y,  y',  z*,  t°)w(x',  y',  z*,  t°)  dx'  dy' 

G(x,  x’ ,  y,  y',  z*,  f°)w(x',  y',  0,  t°)  dx '  dy ' 

(3.5) 


5d>(x,  y,  z',  f°)  , 

-  dz  , 


dz' 


where  the  Green  function  G  and  the  vertical  gradient  of 
the  poloidal  potential  ddVdz  are  given  by 

G(x,  x',  y,  y',  z,  t°) 

M 

=  2  /A>1d)„,(x,  y,  z,  f°)4>„,(x',  y',  z,  f°),  (3.6) 

«j=  1 


d<J>(x,  y,  z,  f°) 
dz 


M 


2  f°)d)m(x,  y,  z). 


(3.7) 


Equation  (3.5)  shows  that  to  determine  w,  the  vertical 
gradient  of  poloidal  potential  should  be  known  on  all 
horizons  from  surface  to  z*  or  from  z*  to  bottom. 

The  same  approach  can  be  applied  for  the  vertical 
vorticity: 


G(x,  x',  y,  y',  z,  f°)^(x',  y',  z,  t°)  dx'  dy' 

=  V(x,  y,  z,  t°),  (3.8) 


where 

G(x,  x',  y,  y',  z,  t°) 

K 

=  2  V^x,  y,  z,  t°)xVk(x\  y',  z,  f°),  and  (3.9) 

A:=  1 


d'Xx,  y,  z,  t°)  =  2  ak(.z,  t°)Vk(x,  y,  z,  t°). 

Jc=l 


(3.10) 


The  series  (3.6),  (3.7),  (3.9),  and  (3.10)  converge  rap¬ 
idly. 


c.  Several  steps  in  reconstruction 

The  reconstruction  is  generally  divided  into  several 
steps:  (a)  parameterizing  k  for  the  unknown  open 
boundary  condition,  (b)  reducing  the  optimization  prob¬ 
lem  (3.2)  to  a  linear  regression  model  for  a  known  open 
boundary  condition  or  a  nonlinear  regression  model  for 
an  unknown  open  boundary  condition,  (c)  determining 
the  optimal  number  of  toroidal  and  poloidal  modes  in 
(3.3),  and  (d)  estimating  the  reconstruction  skill. 

4.  Regularization  technique  for  solving  the 
optimization  problem 

a.  Spectral  expansion  of  k( t) 

We  expand  k(t )  along  the  open  boundary  T'  into 
series: 


S 


k(t)  =  2  csTs(t), 

s=  1 


(4.1) 


where  Ts  (s  =  1,  .  .  . ,  S)  are  the  continuous  polynomials 
of  1,  t,  t2,  .  .  . ,  ts  that  are  preliminarily  orthogonalized 
along  the  open  boundary. 

For  known  (unknown)  boundary  conditions,  a  set  of 
linear  (nonlinear)  regression  equations  is  easily  obtained 
from  (3.2)  to  determine  K  +  M(Q  =  K  +  M  +  S ) 
spectral  coefficients  from  (3.3).  The  nonlinear  regres¬ 
sion  is  discussed  here  only  since  the  nonlinear  regres¬ 
sion  at  each  iteration  uses  the  same  technique  as  the 
linear  regression.  Thus,  it  is  necessary  to  estimate  K  + 
M  spectral  coefficients  in  (3.3)  and  S  coefficients  (c,, 

•  ■  ■ ,  cs)  in  (4.1). 
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b.  Iteration  process  for  nonlinear  regression 

The  classical  iteration  approach  (Eykhoff  1973)  is 
used  to  the  nonlinear  regression  equations.  For  the  ;th 
iteration. 


hIi+n  =  h[i|  +  Shi'1,  (4.2) 

h™  =  (ai/],  .  .  .  ,  a™,  bf,  ...  ,  b$,  cf’,  .  .  . ,  dp)T,  (4.3) 


where  the  increment  <5hm  is  determined  through  mini¬ 
mizing  the  following  functional: 

2 


/(Shi-0  =  ]- 


“°bS  -  M»(hlil)  - 


duk 


hi" 


1 

H — 
2 


'obs  -  u™(hi'J)  - 


ah 

dvKM 


Sh[i] 


ah 


Sh''' 


hi" 


,  (4-4) 


where  the  symbol  (.  .  .)T  indicates  the  transpose  operator. 
The  iteration  process  (4.2)-(4.4)  is  divided  into  nine  steps 
(Fig.  3).  Detailed  description  is  listed  in  appendix  C. 

At  each  iteration,  determining  5hm  is  reduced  to  solv¬ 
ing  a  set  of  ill-posed  linear  algebraic  equations: 

ASh[iI  =  Y,  (4.5) 

where  the  coefficient  matrix  A  is  given  by 


A„  =  B.Bn  +  C.C.. 


Br 


d UKM 

dhr  ’ 


Cr 


9vKm 
dhr  ’ 


r,  q  —  l,  ...,  Q, 

and  the  source  vector  Y  is  represented  by 

Yr  =  [«°bs  -  uKM(h^)]Br  +  [u°bs  -  vKM(h «)]Cr. 

Here,  MATFAB  (1997)  software  is  used  to  compute  the 
basis  functions  {'P*}  and  {<&„,}  and  adopt  the  same 
approach  as  Chu  et  al.  (1997)  to  calculate  Br  and  Cr. 
Obviously,  (4.5)  is  an  ill-posed  system  and  should  be 
solved  using  a  regularization  method. 


c.  Stabilized  least  squares  method 

Numerous  regularization  techniques  exist  such  as  Go¬ 
lub’s  dumping  high-order  singular  values,  Tikhonov’s 
regularization,  iteration  regularization,  and  others  (see 
Engl  et  al.  1996).  Most  existing  regularization  tech¬ 
niques  require  a  priori  knowledge  of  statistical  prop¬ 
erties  of  noise  and/or  the  structure  of  the  solution  of 
(4.5).  Such  a  priori  information  is  usually  not  available 
in  oceanographic  studies.  Therefore,  the  methods  that 
do  not  require  this  a  priori  knowledge  are  useful.  The 
cross-validation  method  (e.g.,  Washba  and  Wendelber- 
ger  1980)  is  very  popular  in  geophysical  studies  (e.g., 
Brankart  and  Brasseur  1996).  However,  it  was  found 
that  the  cross-validation  approach  often  overdetermines 
the  regularization  parameter  (Kugiumtzis  et  al.  1998; 
Ivanov  et  al.  2001a)  and  is  highly  sensitive  to  the  size 
of  observational  samples  used  for  the  reconstruction 
(Mikhalsky  1987;  Ivanov  et  al.  2001a).  Thus,  the  sta¬ 


Model  Inpul 


Fig.  3.  Flowchart  for  determining  the  realization  (hm)  using  the 
iteration  process  (4.2)-(4.4). 


bilized  least  squares  (SFS)  method  is  used  (Ivanov  and 
Margolina  1996;  Ivanov  et  al.  2001a). 

This  method  constructs  a  special  quasi-orthogonal  ro¬ 
tation  matrix  R  and  multiplies  it  by  the  specially  con¬ 
structed  linear  system 

RA'Sh^  =  RY',  (4.6) 

such  that  RA'  — >  Areg  and  RY'  — >  Yreg,  where  Areg  and 
Yreg  are  both  regularized  Q  X  Q  matrix  A  and  Q- di¬ 
mensional  observation  vector  Y.  Appendix  B  shows  the 
procedure  to  construct  the  Q  X  2P  matrix  A'  and  the 
2P-dimensional  observation  vector  Y'.  A  new  trans¬ 
formed  system  (4.6)  should  have  a  much  smaller  noise- 
to-signal  ratio  17,  and  conditional  number  17,  (the  ratio 
of  the  maximum  to  minimum  singular  values  of  the 
system  matrix). 
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The  rotation  matrix  is  constructed  as  the  limit  of  plane 
rotation  matrices  (Givens  matrices): 


J  /  2P  2P  \ 

r  =  I™  n  ww 1  •  ■  ■  n  w  n  w  >  ^ 

J->°°  j=  1  \  P= 3  p= 2  ) 

where  Rf  is  Givens  matrix  at  jth  iteration  applied  to 
q\h  and  pth  lines.  The  rotation  angles  (cp^f)  of  R\'n)  are 
determined  through  the  following  maximization  pro¬ 
cess: 


JAW)  =  2  (W^y-WShJH  -  (R\?Yu-»y  -A  1 


,[0]  = 


yroi  =  y 

p  p’ 


with  ||Sheff|||  called  the  effective  norm. 
This  norm  is  determined  by 


(4.8) 


J  = 


\ll<5heff||2  ) 


— ¥  min, 


(4.9) 


where  5hm  is  a  regularized  solution  of  the  transformed 
system  (4.6).  Iteration  is  used  to  solve  (4.8)  and  (4.9) 
with  an  empirical  requirement: 


1.5  P  <  Q  <  2P. 


See  Ivanov  et  al.  (2001a)  for  detailed  information. 


d.  Properties  of  the  rotation  matrix 

To  construct  the  linear  transformation  R.  we  do  not 
need  a  priori  knowledge  on  statistical  properties  of 
noise.  However,  the  best  results  are  obtained  if  the  exact 
observational  (no  error)  vector  Yexact  is  perpendicular  to 
the  noise  vector  Ynoise;  that  is,  the  inner  product  is  equal 
to  0: 


■  •  I  =11 

exact  noise 

which  indicates  no  interference  between  noise  and  the 
exact  observational  vectors. 

For  a  high  conditional  number,  the  rotation  matrix 
algorithm  (4.6)-(4.9)  has  the  same  accuracy  as  Tikhon¬ 
ov  et  al.’s  (1995)  technique  for  selecting  the  optimal 
regularization  parameter  (Ivanov  et  al.  2001a).  For  a 
small  condition  number  (17,  — >  1),  the  estimate  using 
the  SLS  method  has  the  same  accuracy  as  the  classical 
least  squares  (LES)  solution  (4.4)  when  the  noise-to- 
signal  ratio  is  small  (tj,  — >  0),  and  has  much  better 
accuracy  than  the  LES  estimation  when  the  noise-to- 
signal  ratio  is  large  (tj,  »  0). 

The  accuracy  of  the  approach  is  estimated  by 

cov(Sh)  =  (RTA')-'RTfiR(RTA')-1,  (4.10) 

where  il  is  the  noise  covariance  matrix,  which  is  usually 
unknown,  and  which  makes  (4.10)  less  practical. 


Model  Input 


Fig.  4.  Flowchart  for  determining  the  first  guess  ( 8h|0] )  for  the 
iteration  process  (4.2)-(4.4). 


e.  First  guess  for  the  iteration  process 


It  is  well  known  that  the  iteration  process  (4.2)-(4.4) 
(shown  in  Fig.  3)  may  be  divergent  if  the  first  guess 
Shm  is  far  from  reality.  Therefore,  to  have  a  good  es¬ 
timate  of  (5h|n|  becomes  important.  A  boundary  exten¬ 
sion  method  is  proposed  to  determine  the  initial  guess 
values  of  the  spectral  coefficients.  Let  the  domain  (12,) 
have  an  open  boundary  F[.  The  domain  12,  is  extended 
into  12,  (12,  C  12,)  with  an  open  boundary  T'2  (Fig.  1) 
such  that 


dn 


r;  =  o. 


12 


d<F 

dZ 


=  0. 


(4.11) 


The  circulation  is  reconstructed  in  12,  using  (4.5)  from 
appropriate  linear  regression  model  at  first.  Then,  with 
the  known  velocity  field  at  T[,  we  determine  k  using 
(2.9)  and  obtain  <5hl()|  for  the  domain  12,.  Figure  4  shows 
the  flowchart  for  determining  bh10  .  Detailed  description 
of  the  process  is  listed  in  appendix  D. 
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f.  A  priori  limitation  in  the  regularization  procedure 

The  construction  of  quasi-orthogonal  matrix  R  uses 
an  a  priori  limit  formulated  as  an  inequality.  Strakhov 
(1991)  pointed  out  that  in  the  maximizing  process  (4.8), 
the  observational  vector  is  projected  onto  a  hyperplane 
such  that  the  Cauchy-Schwartz-Buniakowsky  inequal¬ 
ity  is  satisfied: 

i  A'prA'rqm%  -  ( Y'p )2  >  0,  (4.12) 

q=  i 

where  A  and  Y  are  the  system  matrix  and  observational 
vector  after  the  transformation  R  is  applied. 

Such  a  procedure  reduces  the  noise-to-signal  ratio  in 
the  right-hand  side  of  the  transformed  algebraic  equa¬ 
tions  (Eykhnoff  1973;  Strakhov  1991;  Ivanov  et  al. 
2001a)  and  is  called  the  “filtration  of  linear  algebraic 
equations”  (Strakhov  1991).  Notice  that  the  inequality 

(4.12)  is  automatically  satisfied  when  the  observation 
is  not  distorted  by  noise. 

Another  way  to  introduce  an  a  priori  limitation  is  to 
obtain  a  highly  smoothed  solution  8hsmoth  at  first.  Such 
a  solution  is  stable  to  perturbations  with  different  scales 
and  its  main  peculiarities  are  already  lost.  Then,  the 
regularized  solution  8 h  is  found  in  the  neighborhood  of 
this  smoothed  solution  8lismoth  by  adding  a  term, 
y||8/!smoth  —  8h\\P,  to  the  cost  function.  However,  it  has 
not  been  used  in  the  current  study. 


g.  Vertical  velocity 

After  reconstructing  the  horizontal  velocity  field  at 
various  depths,  the  integral  equation  (3.5)  is  solved  nu¬ 
merically  to  obtain  vertical  velocity  by  discretizing  into 
a  set  of  algebraic  equations  in  the  ordinary  way  (Tik¬ 
honov  et  al.  1995): 


method  cannot  be  used  to  regularize  the  ill-posed  system 

(4.13). 

Tikhonov  et  al.  (1995)  change  (4.13)  into  an  opti¬ 
mization  process: 

j|Gw  —  z||2  +  S||w||2  — ¥  min,  (4.14) 

where  8  is  the  regularization  parameter.  Here,  the  min¬ 
imum  sensitivity  of  the  solution  to  8  is  used  to  determine 
its  value  without  the  knowledge  of  the  low-order  noise 
statistics.  Baglai  (1986)  pointed  out  that  this  method  is 
equivalent  to  classical  Tikhonov’s  method  in  determin¬ 
ing  8. 


h.  Summation  of  the  Fourier  series 


Usually,  the  calculated  Fourier  coefficients  [ak(z,  t)\, 
{bm(z,  f)l  contain  errors.  If  the  errors  are  large,  the 
summation  of  the  Fourier  series  in  (3.6),  (3.7),  (3.9), 
and  (3.10)  needs  to  be  regularized  [taking  (3.7)  as  an 
example].  Such  a  regularization  is  achieved  through 
multiplying  each  term  of  the  series  by  a  weight  factor 
(Strahkov  and  Valyashko  1981): 


5$ 

dz 


2 


/r2(z)  +  SA2(z) 


bm{z,  t)<VJx,  y,  z), 


(4.15) 


where  A„,  is  chosen  as 


A: 


I  bj 

max(IZ?J) 


Such  a  process  strongly  smooths  out  the  low-energy 
modes  and  keeps  the  high-energy  modes  in  (4.15).  Here¬ 
in,  also,  8  is  estimated  from  sensitivity  ideas. 


5.  Estimation  of  reconstruction  skill  and  choice  of 
model  parameters  (K,  M,  S ) 


Gw  =  z, 


(4.13)  a.  Vapnik-Chervonkis  cost  function 


where  G  is  a  coefficient  matrix,  and  w  and  z  are  vectors. 
Their  dimensions  are  determined  by  the  discretization 
of  (3.5).  Since  G  is  usually  a  square  matrix,  the  SFS 


Theoretically,  the  reconstruction  skill  should  be  es¬ 
timated  through  the  statistical  cost  function  ( J ),  deter¬ 
mined  as 


<7)  =  2 


(m  —  uKM)2ffu,  v\x,  y)f2(x,  y)  du  dv  dx  dy 


+  (v  —  vKMff\  (u,  v\x,  y)f2(x,  y)  du  dv  dx  dy 


(5.1) 


where  /, (u,  v  \  x,  y)  and  f  2(x,  y)  are  both  the  proba¬ 
bility  density  functions  for  the  velocity  and  station  dis¬ 
position  inside  the  domain,  respectively;  the  angle 
brackets  (( ))  denote  the  ensemble  average  over  real¬ 
izations  of  u  and  v. 


For  a  limited  number  of  data  ( P  finite),  the  ensemble 
average  is  replaced  by  the  empirical  average  over  all  ob¬ 
servational  points.  This  leads  to  the  empirical  cost  function 


7emP  =  J(a\,  ■  ■■,  aK,  bi, ... ,  bM,  k,  P).  (5.2) 
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For  a  large  number  of  data  (P  — >  °°),  the  coefficients 
{dp,  ... ,  aK,  bu  ... ,  bM,  k]  obtained  by  minimizing  Jemp 
coincide  with  the  coefficients  {af,  .  .  .  ,  af,  b f,  .  .  . ,  b%, 
k*  }  obtained  by  minimizing  (J).  However,  this  is  not  true 
for  the  limited  data. 

Vapnik  (1983)  suggested  to  estimate  the  difference 
between  J  and  (J)  by  the  probabilistic  approach.  He 
found  that  the  probability  of  the  maximum  difference 
between  the  empirical  and  ensemble  cost  functions  be¬ 
ing  larger  than  a  given  tolerance  /j,  should  be  less  than 
some  level  g,  which  tends  to  zero  as  the  observational 
sample  size  tends  to  infinity: 

Probj  sup  \(J(K,  M,  S ))  —  Jemp(K,  M,  5)1  >  /a  l 

I  K,M,S 

^  g(P,  p)  (5.3) 

lim  g(P,  /a )  =  0,  (5.4) 

P—>CC 

where  “Prob”  denotes  the  probability.  From  (5.3)  we 
have 

Jemp(K,  M,  5)  -  i±  <  (J(K,  M,  5)> 

<  Jemp(K,  M,  5)  +  n,  (5.5) 

which  leads  to  a  formula  for  so-called  Vapnik-Cher- 
vonkis  cost  function  (VCCF): 


1 

sup  ~[(uobs  -  uKmy 

ai,...,aK,b\,...,bM,K,x,y  ^ 

+  (”°bs  -  vKJ2]  s  L  (5.8) 

Then,  the  optimal  parameter  Mopt  should  be  obtained 
through  the  minimization  process: 

^optMopt  _  min(^s:0p10’  Jk0 p,l>  J Kopt2 ’  •  •  •  >  J kpp1m)-  (5-9) 

The  reconstructed  field  {u(Kopt,  Mopt),  v(Kopt,  Mopt)} 
is  the  solution  that  guarantees  the  minimum  difference 
between  (J)  and  /emp  for  a  given  observational  sample 
size.  It  is  noticed  that  the  less  the  difference  between 
(J)  and  /emp,  the  closer  the  truncated  mode  {uKM,  vKM] 
is  to  reality. 

From  the  physical  point  of  view,  the  optimal  param¬ 
eters  are  found  through  the  cost  function  (5.9),  which 
is  a  trade-off  between  the  likelihood  of  a  model  (i.e., 
its  ability  to  reproduce  data)  and  a  penalty  for  making 
the  model  too  complex  for  the  data. 

c.  Root-mean-square  error 

When  the  exact  solution  is  known,  we  used  the  root- 
mean-square  error  (rmse)  inside  the  domain: 


Jkms  =  Jcmp(K,  M,  S )  +  /A.  (5.6) 

The  tolerance  /a  depends  on  the  chosen  numerical 
model  and  is  determined  by  two  criteria:  (a)  minimi¬ 
zation  of  the  probability  (5.3),  and  (b)  uniform  con¬ 
vergence  of  the  coefficient  set 

{«i.  .  .  .  ,  dK,  bu  .  .  . ,  bM,  i<} 

to  {af,  .  .  . ,  a%,  bf,  .  .  .  ,  b%,  k*}  as  P  — > 


b.  Optimal  spectral  truncation 

Vapnik’s  (1983)  method  is  used  to  determine  the  op¬ 
timal  values  for  three  model  parameters  (Kopt,  Mopt,  Sopt). 
For  simplicity  without  loss  of  generality,  we  illustrate 
the  approach  with  a  given  5opt.  Let  the  family  of  models 

(uK0,  vK0),  (uK ,,  vKl),  (uK 2,  VK2),  •  •  .  ,  ( UKM ,  vKM ) 


Xv  = 


and  along  the  open  boundary: 


Xn  = 


(V„  ~  V,,)2 

(V2„) 


dr' 


(5.10) 


to  determine  the  reconstruction  skill.  Here,  wex,  u.x,  V,„ 
and  V T  represent  the  exact  velocity  components,  (')  de¬ 
notes  the  reconstructed  field,  and  (.  .  .)  is  the  averaging 
operator  over  the  domain  or  along  the  open  boundary. 
To  estimate  spatial  structure  of  reconstruction  error,  the 
residual  velocity  defined  by  Lipphardt  et  al.  (2000), 


be  used  in  the  reconstruction  process.  The  velocity  com¬ 
ponent  (uK0,  vK0)  does  not  contain  contribution  from  the 
poloidal  potential  (i.e.,  b„  =  0).  The  optimal  parameter 
K  t  is  determined  through  minimization  of  the  func¬ 
tional 


„  \K 

^ Km  ^emp  ~  'A,  j  p 


m  =  0,  .  . . ,  M, 


(5.7) 


ures  =  uex  -  nKM ,  (5.12) 

is  used.  It  is  noticed  that  such  a  skill  is  on  the  base  of 
the  known  exact  velocity  uex  and  ( V n ,  VT)  and  has  little 
practical  significance.  However,  it  is  applicable  for  ver¬ 
ification  in  identical  twin  experiments. 

6.  Computational  cost 

Flow  reconstruction  is  performed  by  solving  a  set  of 
(K  +  M)I  linear  algebraic  equations  (3.3)  for  the  zero- 
order  approximation  case  and  by  solving  a  set  of  ( K  + 


where  g  =  0.95  and 
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M  +  S){K  +  M  +  S)IIt  +  (K  +  M  +  5)//,  linear 
algebraic  equations  (4.1)-(4.5)  for  iteration  process  for 
nonlinear  regression.  Here,  I  and  /  are  the  numbers  of 
iteration  to  determining  h[i]  and  (ATopt,  Mopt),  respective¬ 
ly.  Thus,  computational  cost  depends  on  the  selection 
of  the  optimal  modes  (the  basis  functions)  used  in  the 
spectral  decomposition.  In  general,  the  number  of  the 
basis  functions  determines  the  smoothness  of  the  re¬ 
constructed  circulation  (Courant  and  Hilbert  1966). 

In  additional  to  the  three  examples  presented  in  Part 
II  of  this  paper,  the  spectral  decomposition  (2. 1 1  )-(2. 17) 
was  performed  for  the  Black  Sea  with  horizontal  res¬ 
olution  of  18  km  (Eremeev  et  al.  1992a, b);  for  the  Gulf 
of  Mexico  with  horizontal  resolution  of  9  km  (Chu  et 
al.  2003);  for  the  Monterey  Bay  with  horizontal  reso¬ 
lution  of  1  km  (Ivanov  and  Melnichenko  2002);  and  for 
Louisiana-Texas  shelf  with  horizontal  resolution  of  10 
km  (Ivanov  et  al.  2001b).  In  these  calculations,  the  op¬ 
timal  number  of  basis  functions  are 

Kopt  =  40,  Mopt  =  30,  S  =  0, 

for  the  Black  Sea  reconstruction  [closed  basin,  with  the 
boundary  conditions  (2.11),  (2.12)]; 

Kopt  =  120,  Mopt  =  20,  Sop{  =  5, 

for  the  Gulf  of  Mexico  reconstruction  [semiclosed  basin, 
with  the  open  boundary  conditions  (2.10),  (2.13)]; 

Kopt  <  70,  Mopt  <  20, 

for  the  Monterey  Bay  reconstruction  [with  the  open 
boundary  conditions  (2.10),  (2.13),  and  (2.18)];  and 

Kopt  <  40,  Mopt  <  10,  Sopt  =  3, 

for  the  Louisiana-Texas  shelf  reconstruction  [open  shelf 
area,  with  the  simplified  open  boundary  conditions 
(2.18)].  If  the  zero-order  approximation  (first  guess)  is 
successfully  chosen,  the  numbers  of  iteration,  I  and  /, , 
usually  do  not  exceed  7.  If  they  are  larger  than  7,  the 
simplified  open  boundary  conditions  (2.18)  is  recom¬ 
mended. 

In  reconstructing  the  Gulf  of  Mexico  circulation  on 
a  Pentium-Ill  (850  MHz)  processor,  300  toroidal  and 
poloidal  basis  functions  are  determined  (40  CPU  min) 
from  a  15  400  X  15  400  matrix  corresponding  to  2D 
Laplacian  with  the  Gulf  of  Mexico  geometry.  A  set  of 
105  linear  algebraic  equations  is  solved  at  each  iteration 
in  (4.2)-(4.5)  using  a  special  technique  for  super  large 
linear  algebraic  system  developed  by  Strakhov  and 
Strakhov  (1999).  Such  calculation  consumes  about  5.6 
CPU  h.  The  relatively  high  computational  cost  is  caused 
by  the  use  of  large  data  with  the  total  open  boundary 
conditions  (2.10)  and  (2.13).  If  the  simple  open  bound¬ 
ary  condition  (2.10)  is  used,  the  computational  cost  is 
greatly  reduced. 

7.  Conclusions 

First,  a  multistep  scheme  is  developed  to  reconstruct 
velocity  from  sparse  and  noisy  data  in  an  open  domain: 


(a)  a  boundary  extension  method  to  determine  normal 
and  tangential  velocities  at  an  open  boundary,  (b)  es¬ 
tablishment  of  homogeneous  open  boundary  conditions 
for  the  two  potentials  with  a  spatially  varying  coefficient 
k(t),  (c)  spectral  expansion  of  k(t),  (d)  determination 
of  basis  functions  for  the  two  potentials  for  the  spectral 
decomposition  using  homogeneous  boundary  condi¬ 
tions,  and  (e)  determination  of  coefficients  in  the  spec¬ 
tral  decomposition  of  velocity  and  k(t)  using  linear  or 
nonlinear  regressions. 

Second,  the  homogeneous  boundary  conditions  of  (Mf, 
d>)  at  both  rigid  and  open  boundary  segments  make  it 
possible  to  obtain  basis  functions  for  an  open  domain. 
The  basis  functions  are  the  eigenfunctions  of  the  La¬ 
placian  operator  with  homogeneous  boundary  condi¬ 
tions  and  depend  on  the  spectrally  varying  parameter  k 
at  the  open  boundary. 

Third,  the  spectra  of  the  two  horizontal  velocity  com¬ 
ponents  and  k  are  truncated  at  K,  M,  and  S.  The  optimal 
model  parameters  ( Kopt ,  Mopt,  .S-op, )  are  determined 
through  a  modified  cost  function,  which  is  constructed 
on  the  basis  of  model  capability  and  data  reproduction 
complexity  (penalty).  This  cost  function  is  also  used  to 
verify  the  model  reconstruction  skill  from  sparse  and 
noisy  data. 

Fourth,  the  spectral  coefficients  for  horizontal  veloc¬ 
ity  and  k  are  determined  simultaneously  using  the  sta¬ 
bilized  least  squares  (SLS)  method.  This  method  does 
not  require  a  priori  knowledge  about  noise  and  is  robust 
to  the  size  of  observational  samples  used  for  the  re¬ 
construction. 

Fifth,  after  reconstructing  the  horizontal  velocity  field 
at  various  depths,  the  vertical  velocity  may  be  recon¬ 
structed  through  solving  the  integral  equation  (3.5)  nu¬ 
merically.  Since  the  coefficient  matrix  is  square,  the 
minimum  sensitivity  of  solution  is  used  to  determine 
the  regularization  parameter  and  then  use  Tikhonov’s 
approach  to  reconstruct  the  vertical  velocity. 
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APPENDIX  A 

Two-Scalar  Potential  Representation  of 
Three-Dimensional  Incompressible  Flow 

Consider  spherical  coordinates  with  r  representing  the 
radius  vector.  Using  the  Helmholtz  theorem  (Morse  and 
Feshbach  1953),  the  horizontal  velocity  at  a  given  radius 
I  r  I  =  r  is  written  as 


rA2  +  VA3  =  V  X  V  X  (r<J>).  (A.8) 

Let  A,  =  ■'P,  then  (A. 3)  becomes 

u  =  V  X  (DP)  +  V  X  V  X  (r<I>),  (A. 9) 

which  is  the  two-scalar  representation  for  a  three-di¬ 
mensional  incompressible  field.  For  more  details,  see 
Moffat  (1978)  and  Zeldovich  et  al.  (1985).  Obviously, 
the  two  potentials  P'P,  <T>)  defined  here  are  not  the  same 
as  the  streamfunction  (A , )  and  velocity  potential  (A3) 
of  a  two-dimensional  flow  on  a  spherical  surface. 


u„  =  r  X  V„A,  +  V„A3,  (A.l) 

where  VH  is  the  two-dimensional  nabla  operator. 

Let  us  define  a  new  function  A,  such  that  the  radial 
velocity  component  is  written  by 

dA3 

iir  =  rA 2  H - -.  (A. 2) 

dr 


Combination  of  (A.l)  and  (A. 2)  gives  a  well-known 
form  for  a  three-dimensional  vector  field  (Dubrovin  et 
al.  1992): 


u  =  r  X  VA,  +  rA2  +  VA3.  (A. 3) 

To  reduce  the  number  of  scalar  functions  in  (A.  3) 
from  3  to  2,  two  fundamental  properties  of  oceanic  flows 
are  accounted  for  1)  existence  of  the  vertical  axis 
marked  out  by  the  earth’s  rotation  for  large-,  meso-  and 
submesoscales;  and  2)  incompressibility. 

The  first  property  of  ocean  flows  allows  to  construct 
the  unique  invariant  transformation  for  (A.l)  and  (A.2): 


A,  =  A,  -  Mr), 


a2  =  a2  +  !^, 

r  dr 


A3  =  A3  -  f2(r),  (A. 4) 

where  /,  and  f2  are  two  arbitrary  functions  depending 
on  r  only.  In  general,  the  three-dimensional  flow  is  de¬ 
composed  by 

u  =  r  X  VA,  +  rA2  +  VA3.  (A.5) 


APPENDIX  B 


Structure  of  A '  and  Y'  in  (4.5) 

The  matrix  A'  and  the  vector  Y'  in  (4.6)  are  given  by 


A' 


du , 


d/z. 


du. 


d/z, 
du ... 


8hl 


du. 


d/z. 


p=  i 


P=p 


p=  i 


P=p 


du. 


d  hr 


du. 


d/z. 


du. 


d  hr. 


dv» 


d/z. 


p=  1 


p=p 


p=  1 


p=p 


Y'  = 


W “bs  -  ^M(h[,1)W^ 

«ubs  —  UKMW])\p=P 
W°bs  -  uKM(h[,1)lp=i 

V°bs  ~  VKJhMP-Pj 


(B.l) 


Two  velocity  components, 

utor  =  r  X  VA,  =  V  X  (rA,),  and 

up°i  =  M2  +  VA3,  (A. 6) 

are  also  invariant  in  the  transform  (A.  3).  The  two  com¬ 
ponents  utor  and  upo1  are  called  the  toroidal  and  poloidal 
velocities  (Moffat  1978),  respectively.  Substitution  of 
(A. 3)  into  the  incompressibility  condition 

V  u  =  0 

leads  to  the  following  relationship: 

AA3  =  —  V(rA2).  (A. 7) 

According  to  (A. 6),  the  poloidal  potential  <f>  of  the  field 
can  be  defined  by 


After  multiplying  R  to  the  system  (4.5)  we  have 


^11 

a1Q 

©  ^ 

b,  ' 

© 

aQl 

^QQ 

© 

;  Y  = 

5  reg 

bQ 

© 

0  0 

© 

0 

VC2P1 

C2P2P  J 

M2P2P, 

(B.2) 

Here,  the  last  rows  in  A,,,.,  and  Yreg  are  excluded  from 
the  analysis  because  they  are  not  filtered.  The  trans¬ 
formed  system  (4.6)  is  well  posed  and  can  be  solved 
using  a  typical  linear  algebraic  method  such  as  the  Gauss 
procedure. 
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APPENDIX  C 

Determination  of  the  Realization  (Hl/1) 

This  process  starts  from  the  first-guess  optimal  ve¬ 
locity  OgplMopt,  <i)iAV)  in  the  domain  11,  (T  U  T[).  There 
are  nine  steps  (for  iteration  ;)  to  determine  the  reali¬ 
zation  (hm)  (Fig.  3). 

1)  Calculate  k1'1  (i  =  1,  2,  .  .  .)  along  the  open  boundary 
TJ  for  the  domain  (1  (T  U  TJ). 

2)  Determine  the  basis  functions  {'TY1,  \FY]>  .  .  . , 

and  { (fi1,'1,  ,  d>0]}  for  the  domain  D,  with 

K[il. 

3)  Calculate  the  velocity  (ul^]M,  vl£h)  for  the  pair  of  (K, 
M).  Determine  the  optimal  pair  (Xopt,  Mopt)  through 
minimizing  the  Vapnik-Chervonkis  cost  function 
(5.6).  Then,  obtain  the  optimal  adjustment  6h^|  iMp 
using  (4.6). 

4)  Update  the  vector  hm. 

5)  Compute  the  optimal  velocity  (M^pt„opi,  u£ptM  ). 

6)  Compare  ||Sh|,|||  with  the  tolerance  level  e.  If  ||Shm|| 

<  e,  the  iteration  stops,  and  (w^plMopl,  )  is  the 

final  reconstructed  velocity  field.  If  ||Sh[il||  >  e,  the 
next  step  iteration  starts. 

APPENDIX  D 

Determination  of  the  First  Guess  (5H|01) 

This  process  starts  from  noisy  observational  velocity 
data  in  the  expanded  domain  12  2  (T  U  T2)  with  the 
assumption  Km  =  0.  There  are  five  steps  to  determine 
the  first  guess  (<5h101)  (Fig.  4). 

1)  Calculate  the  two  sets  of  the  basis  functions  {"'F1,01, 

4,i201,  .  .  . ,  ''pm}  and  {€>m  .  .  .  ,  <&$}  using  the 

Lanchoz-Arnoldi  method  from  MATLAB  6.0. 

2)  Calculate  the  velocity  vlgh)  for  the  pair  of  (K, 
M)  using  (3.3).  Determine  the  optimal  pair  (Kopt, 
Mopt)  through  minimizing  the  Vapnik-Chervonkis 
cost  function  (5.6).  Detailed  explanation  of  this  pro¬ 
cedure  can  be  found  in  section  5b. 

3)  Compute  the  first-guess  optimal  velocity  (w^J  M  t, 

using  (3.3). 
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